Tarea 4: Bayesian Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

#Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
library(purrr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
# Análisis bayesiano
library(rethinking)
## Loading required package: rstan
## Loading required package: StanHeaders
## 
## rstan version 2.32.6 (Stan version 2.32.2)
## 
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## 
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Loading required package: parallel
## Loading required package: dagitty
## rethinking (Version 2.01)
## 
## Attaching package: 'rethinking'
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")
## Warning: package 'rethinking' is in use and will not be installed

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice el siguiente chunk:

install.packages(c("mvtnorm","loo","coda", "shape"), repos="https://cloud.r-project.org/",dependencies=TRUE)
## Installing packages into 'C:/Users/juanm/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)
## package 'mvtnorm' successfully unpacked and MD5 sums checked
## Warning: cannot remove prior installation of package 'mvtnorm'
## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying
## C:\Users\juanm\AppData\Local\R\win-library\4.4\00LOCK\mvtnorm\libs\x64\mvtnorm.dll
## to C:\Users\juanm\AppData\Local\R\win-library\4.4\mvtnorm\libs\x64\mvtnorm.dll:
## Permission denied
## Warning: restored 'mvtnorm'
## package 'loo' successfully unpacked and MD5 sums checked
## package 'coda' successfully unpacked and MD5 sums checked
## package 'shape' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\juanm\AppData\Local\Temp\RtmpCyopNW\downloaded_packages
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')
## Warning: package 'rethinking' is in use and will not be installed

Pregunta 1:

Un conjunto de carteros aburridos de las mordidas de perros ha decidido realizar un catastro de mordidas recibidas por los empleados de su empresa en un periodo de 8 meses, planeando en base a estos datos realizar inferencia bayesiana. Los datos de las mordidas estas datos por el dataset no+mordidas.csv, en donde cada fila representa las mordidas recibidas por diferentes carteros y las columnas señalan si fue mordido o no el cartero en los meses de estudio (si fue mordido es señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar que un cartero no puede ser mordido mas de una vez al mes, ya que el damnificado recibe licencia de todo un mes tras la mordida, reincorporándose el siguiente mes al trabajo.

dataMordidas <- read.csv("no_mordidas.csv", header = TRUE)
cols <- c("bites_month_1", "bites_month_2", "bites_month_3",
          "bites_month_4", "bites_month_5", "bites_month_6",
          "bites_month_7", "bites_month_8")

# Prior parameters
grid_len <- 250
p_grid <- seq(from=0, to=1, length.out = grid_len)
prior <- dnorm(p_grid, mean = 0.6, sd = 0.05)

posteriors <- c()
counted.months <- c()
p_grids <- rep(p_grid, 4)

for (n in c(1, 2, 4, 8)) {
  bites <- rowSums(as.matrix(dataMordidas[, 1:n]))

  likelihood <- dbinom(bites, size = n, prob = p_grid)
  
  unstd.posterior <- prior * likelihood
  posterior <- unstd.posterior / sum(unstd.posterior)
  
  posteriors <- c(posteriors, posterior)
  counted.months <- c(counted.months, rep(n, grid_len))
}

post.df <- data.frame("grid" = p_grids, "posterior" = posteriors, "months" = counted.months)

# Visualize
ggplot(post.df) +
  geom_line(aes(x = grid, y = posteriors)) +
  facet_grid(months ~ .) +
  labs(x = "Probabilidad de Mordida", y = "Posterior", title = "Distribución Posterior por Número de Meses")

Respuesta:

Es difícil de interpretar por como está graficado. Sin embargo, podemos ver que el comportamiento del posterior a medida que avanzan los meses se ve fuertemente influenciado por el prior. Vemos que muestra un comportamiento variable con una deviación cercana a la dada por el carteror (0.05) y asociado a una probabilidad de mordida del 0.6.

cols <- c("bites_month_1", "bites_month_2", "bites_month_3",
          "bites_month_4", "bites_month_5", "bites_month_6",
          "bites_month_7", "bites_month_8")

# Prior uniforme
grid_len <- 250
p_grid <- seq(from=0, to=1, length.out = grid_len)
prior <- rep(1, grid_len)

posteriors <- c()
counted.months <- c()
p_grids <- rep(p_grid, 4)

for (n in c(1, 2, 4, 8)) {
  bites <- rowSums(as.matrix(dataMordidas[, 1:n]))

  likelihood <- dbinom(bites, size = n, prob = p_grid)
  
  unstd.posterior <- prior * likelihood
  posterior <- unstd.posterior / sum(unstd.posterior)
  
  posteriors <- c(posteriors, posterior)
}

post.df$posterior2 <- posteriors

# Graficar los resultados
ggplot(post.df) +
  geom_line(aes(x = grid, y = posteriors)) +
  facet_grid(months ~ .) +
  labs(x = "Probabilidad de Mordida", y = "Posterior", title = "Distribución Posterior por Número de Meses")

Respuesta:

Vemos que el comportamiento anterior no se ve reflejado del todo. Aquí, la distribución del posterior parece mantenerse con el avance de los meses. Esto podría indicar que es importante la elección del prior, ya que, en este caso al menos, la elección de este cambia los resultados del análisis Bayesiano.

for (n in c(1, 2, 4, 8)) {
  
  W <- sum(dataMordidas[, 1:n])
  L <- n * 250 - W
  
  # Adjust model
  res <- quap(
    alist(
      W ~ dbinom(L, p),  # Likelihood
      p ~ dunif(0, 1)    # Prior 
    ),
    data = list(W = W, L = L)
  )
  
  df <- precis(res)
  
  # Graph the posterior curves
  curve(dnorm(x, df$mean, df$sd), lty=2, add=FALSE, 
        xlab = "Probabilidad de mordida", ylab = "Densidad",
        main = paste("Distribución posterior - Meses:", n))
}

Respuesta:

Los resultados son más interpretables, ya que se observan directamente las densidades de las probabilidades. En este caso, Vemos que efectivamente suben las probabilidades respecto a su valor inicial (0) pero además tenemos un rango más cercano en el que estimar el cual es entre 0.6 y 0.7, dado que la distribución del posterior parece ser una gaussiana en ese rango. Vemos además que con el avance de los meses, el centro de la distribución se va a acercando a 0.6, que corresponde al centro del prior.

¿Cómo puede interpretar los resultados?

grid_len <- 1000
p_grid <- seq(0, 1, length.out = grid_len)

prior <- dnorm(p_grid, mean = 0.6, sd = 0.05)

bites <- rowSums(as.matrix(dataMordidas[, 1:8]))

likelihood <- dbinom(bites, size = 8, prob = p_grid)

unstd.posterior <- prior * likelihood

posterior <- unstd.posterior / sum(unstd.posterior)

samples <- sample(p_grid, prob = posterior, size = 1e4, replace = TRUE)

dens(samples)

# Intervalos:
proporcion_0_035 <- mean(samples >= 0 & samples <= 0.35)
proporcion_035_045 <- mean(samples > 0.35 & samples <= 0.45)
proporcion_mayor_045 <- mean(samples > 0.45)

# Print proportions
cat("Proporción en el rango [0, 0.35]:", proporcion_0_035, "\n")
## Proporción en el rango [0, 0.35]: 0
cat("Proporción en el rango (0.35, 0.45]:", proporcion_035_045, "\n")
## Proporción en el rango (0.35, 0.45]: 0.0022
cat("Proporción en el rango (>0.45):", proporcion_mayor_045, "\n")
## Proporción en el rango (>0.45): 0.9978

Respuesta:

Los resultados parecen indicar una proporción mayor para la posterior cercana a 0.6. La proporción de los defined boundaries indica que la mayor parte de estas posteriores están en el rango mayor a 0.45. Esto nos indica que la probabilidad de mordida es mayor a este valor y es cercana a 0.6.

-[] Calcule un intervalo de credibilidad al \(50%\), \(75%\) y \(95%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?

# Function to calculate credible interval for given prob
cred_interval <- function(samples, prob) {
  quantiles <- quantile(samples, probs = c((1 - prob) / 2, 1 - (1 - prob) / 2))
  return(quantiles)
}

# Calculate intervals
cred_50 <- cred_interval(samples, 0.50)
cred_75 <- cred_interval(samples, 0.75)
cred_95 <- cred_interval(samples, 0.95)

Intervalo \(50\%\):

cat("Intervalo de Credibilidad 50%: [", cred_50[1], ", ", cred_50[2], "]\n")
## Intervalo de Credibilidad 50%: [ 0.5585586 ,  0.6236236 ]

Intervalo \(75\%\)

cat("Intervalo de Credibilidad 75%: [", cred_75[1], ", ", cred_75[2], "]\n")
## Intervalo de Credibilidad 75%: [ 0.5365365 ,  0.6476476 ]

Intervalo \(95\%\)

cat("Intervalo de Credibilidad 95%: [", cred_95[1], ", ", cred_95[2], "]\n")
## Intervalo de Credibilidad 95%: [ 0.4964965 ,  0.6846847 ]

Respuesta:

Los resultados nos indican con importante certeza que la probabilidad se encuentra al menos en el rango entre 0.4 y 0.6. Un problema que pueden tener los intervalos de credibilidad es que estos dependen fuertemente de la elección del modelo, de los datos y de la elección del prior. En este caso, los intervalos de credibilidad se están viendo afectados por la elección del prior, ya que este es centrado en 0.6. Claramente los intervalos tienden a incluir este valor.

# Calculate intervals with HPDI function from rethinking
hdpi_50 <- HPDI(samples, prob = 0.50)
hdpi_75 <- HPDI(samples, prob = 0.75)
hdpi_95 <- HPDI(samples, prob = 0.95)

Intervalo \(50\%\):

cat("Intervalo HDPI 50%: [", hdpi_50[1], ", ", hdpi_50[2], "]\n")
## Intervalo HDPI 50%: [ 0.5645646 ,  0.6266266 ]

Intervalo \(75\%\)

cat("Intervalo HDPI 75%: [", hdpi_75[1], ", ", hdpi_75[2], "]\n")
## Intervalo HDPI 75%: [ 0.5375375 ,  0.6476476 ]

Intervalo \(95\%\)

cat("Intervalo HDPI 95%: [", hdpi_95[1], ", ", hdpi_95[2], "]\n")
## Intervalo HDPI 95%: [ 0.4994995 ,  0.6866867 ]

Respuesta:

Los intervalos son un poco más estrechos, pero sigue sin ser una diferencia demasiado sustantiva. La diferencia entre los intervalos es que los HDPI no asignan igual densidad de probabilidad en las colas.

bites <- rowSums(as.matrix(dataMordidas[, 1:8]))
n <- nrow(dataMordidas)
ratio.real <- sum(bites) / (n * 8)

simulated_bites <- sapply(samples, function(p) {
  rbinom(n, size = 8, prob = p)
})

simulated_ratios <- colSums(simulated_bites) / (n * 8)

ggplot() +
  geom_density(aes(simulated_ratios)) + 
  geom_vline(aes(xintercept = ratio.real), color = "red") +
  labs(x = "Ratio de Carteros Mordidos", y = "Densidad", title = "Distribución de Ratios Simulados vs. Ratio Real")

Respuesta:

Aparentemente, no, el modelo no se ajusta bien a los datos. El ratio real parece estar un poco bajo 0.4 mientras que la distribución de las predicciones está centrada en 0.6 (Que se asocia más al prior).

first_month <- "bites_month_7"
second_month <- "bites_month_8"

# Filter bitten on first month
mordidos_first <- dataMordidas[dataMordidas[[first_month]] == 1, ]

real_recuento <- sum(mordidos_first[[second_month]])

n_mordidos_first <- nrow(mordidos_first)

# Simulate 10000 times
simulated_counts <- sapply(samples, function(p) {
  rbinom(1, size = n_mordidos_first, prob = p)
})

ggplot() +
  geom_density(aes(simulated_counts), alpha = 0.5) +
  geom_vline(aes(xintercept = real_recuento), color = "red", linetype = "dashed") +
  labs(x = "Recuento de Carteros Mordidos", 
       y = "Densidad",
       title = "Distribución de Simulaciones vs. Recuento Real") +
  theme_minimal()

Respuesta:

Ocurre un fenómeno similar al anterior, donde el modelo no se ajustaba bien y centraba la distribución en el prior. Esto probablemente se debe a un error en los cálculos. Es interesante ver sin embargo como el modelo llega a distribuciones similares.

Pregunta 2:

En esta pregunta se trabajara con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.

Usted sabe un dato extra sobre la información, \(\sigma \in [0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\).

# Leer información
data_notas <- read.csv("notas.csv")
data_vector <- as.numeric(data_notas$Notas)

# Función para crear likelihood dado mu y sigma
grid_function <- function(mu, sigma){
   sum(dnorm(data_vector, mean = mu, sd = sigma, log = TRUE))
}

# Valores de la grilla
grid_mu <- seq(1, 7, length.out = 100)
grid_sigma <- seq(0.5, 1.5, length.out = 100)

# Se crea la grilla 2d
data_grid <- expand_grid(grid_mu,grid_sigma)

# Se guarda la likelihod
data_grid$likelihood <- map2(data_grid$grid_mu,data_grid$grid_sigma, grid_function)

# Se transforma el forma de map2 a una columna
data_grid <- unnest(data_grid,cols = c("likelihood"))

# Valores de los priors
prior_mu <- rep(1, length(grid_mu))
# Se resta 0.5 para mantener el argumento de la distribución en [0,1]
prior_sigma <- dbeta(grid_sigma - 0.5, shape1 = 3, shape2 = 6)

# Se crea la grilla 2d de priors
prior <- expand_grid(prior_mu,prior_sigma)

# Se calculan los valores del prior
data_grid$prior <-  map2(prior$prior_mu,prior$prior_sigma, prod)
data_grid <- unnest(data_grid,cols = c("prior"))

# Se calcula el posterior
data_grid$unstd_posterior <-  data_grid$likelihood * data_grid$prior

# Se estandariza el posterior
data_grid$posterior <- data_grid$unstd_posterior/sum(data_grid$unstd_posterior)

# Se ajustan los valores de la posterior para que no sean valores tan pequeños
data_grid$posterior <- (data_grid$posterior - min(data_grid$posterior))/(max(data_grid$posterior)-min(data_grid$posterior))

Respuesta:

Para el parámetro \(\mu\), que representa la media de las notas, se ha seleccionado un prior uniforme en el rango \(([1, 7])\). Al no disponer de información previa que sugiera que ciertos valores de son más probables que otros, es conveniente el uso de esta prior en el intervalo de las notas, de 1 a 7.

Para el parámetro \(\sigma\), correspondiente a la desviación estándar, se ha optado por un prior de distribución beta. La elección se fundamenta en la información adicional proporcionada, que indica que \(\sigma\) se encuentra en el intervalo \(([0.5, 1.5])\) y es más probable que tome valores entre \(([0.5, 1])\) que en \(((1, 1.5])\). Para representar esto, se ha optado los parámetros \(\alpha = 3\) y \(\beta = 8\) para aumentar la probabilidad dentro del rango \(([0.5, 1])\). Graficamente:

plot(grid_sigma, prior_sigma, type = "l", col = "blue",
     main = "Prior para sigma usando Beta(3, 8)",
     xlab = expression(sigma), ylab = "Densidad de probabilidad")

# Añadir una línea horizontal para referencia
abline(h = 0, col = "gray")

# Punto de referencia
valor_x <- mean(data_vector)
valor_y <- 0.8

# Gráfico
punto_comparacion <- tibble(x = valor_x, y = valor_y)

plt <- data_grid %>%
  ggplot(aes(x = grid_mu, y = grid_sigma)) +
  geom_raster(aes(fill = posterior),
    interpolate = T
  )+ 
  geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
  geom_label(
    data = punto_comparacion, aes(x, y),
    label = "Punto Comparación",
    fill = "green",
    color = "black",
    nudge_y = 0, # Este parametro desplaza la caja por el eje y
    nudge_x = 1 # Este parametro desplaza la caja por el eje x
  )+
  scale_fill_viridis_c() +
  labs(
    title = "Posterior para Mean y Standard Deviation",
    x = expression(mu ["Mean"]),
    y = expression(sigma ["Standar Deviation"])
  ) +
  theme(panel.grid = element_blank())

plt

Respuesta:

La región de mayor probabilidad para \(\mu\) se encuentra entre 1 y 3, mientras que para \(\sigma\) la mayor densidad se concentra en el rango de 0.5 a 0.9. Esto sugiere que los datos observados tienden a estar centrados en valores bajos, con una variabilidad moderada. El prior utilizado para \(\sigma\), basado en una distribución Beta(3, 8), favorece valores bajos, lo que se refleja en la forma adecuada del posterior respecto al punto de comparación en el eje de \(\sigma\). Por el contrario, el prior uniforme para \(\mu\) no parece haber sido el mejor, pues según el posterior, la media debería ser cercana al mínimo, cosa que no se corresponde con la media real, que oscila en torno al 6. De todas formas, sin ninguna información adicional, el prior uniforme seguía siendo buena opción, puesto que con algún otro modelo era posible quedar más lejos del valor real y empeorar aún más el resultado.

# Codificamos los datos
x <- 1:length(data_grid$posterior)

# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)

# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]

# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)

# Realizamos las densidades
dens(df)

Respuesta:

El proceso de “sampling from grid approximation” permite obtener muestras representativas del posterior usando la aproximación discreta calculada previamente en la cuadrícula. El posterior combina la verosimilitud de los datos y el prior, pero al usar una aproximación por cuadrícula, se evalúa en puntos discretos. El muestreo transforma esta representación discreta en una muestra continua, proporcionando una visión más precisa de la distribución posterior de los parámetros, lo que permite convertir el posterior discreto en una distribución continua de muestras, lo que facilita el análisis.

En los gráficos generados, se observa en el izquierdo la densidad de las muestras para \(\mu\), que muestra un pico pronunciado alrededor de \(\mu \approx 1.5\), lo que indica que los valores más probables para la media se encuentran en este rango. La distribución es asimétrica y decrece gradualmente a medida que \(\mu\) aumenta, sugiriendo que aunque es posible encontrar medias mayores, estas son menos probables según los datos observados y el prior utilizado, en concordancia con el análisis previo.

En el gráfico derecho, que muestra la densidad para \(\sigma\), se observa un pico alrededor de \(\sigma \approx 0.8\), lo que indica que estos son los valores más probables para la desviación estándar, acorde a la creencia dada. La distribución es aproximadamente simétrica, con una caída rápida después de 1 y un límite inferior en 0.5, reflejando la influencia del prior que restringe \(\sigma\) al intervalo \(([0.5, 1.5]\).

Pregunta 3: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar a aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Para esta pregunta usaremos el dataset Pokemon.csv, que contiene las características de más de 700 personajes.

poke.df <- read.csv("Pokemon.csv")

Respuesta:

Se utiliza un prior de media 100 y desviación estándar 50 para b0, basado en los valores observados en las variables y el conocimiento previo de los rangos en que suelen estar los ataques de los pokemones (entre 0 y 200, siempre positivos). Se asume también que la relación linear de la pendiente es no tan variable, es decir, con una pendiente 1 y variabilidad moderada de 0.5. Finalmente, se elige una desviación estándar conservadora, en el rango 0-200 mencionado anteriormente. Las ecuaciones del modelos son las del código que se encuentran en las alist model.

Además, el código realiza el ajusta mediante la aproximación de Laplace con el método quap, entregando la lista con los parámetros y priors definidos.

# Model variables
data <- poke.df[, c("Sp..Atk", "Attack")]

# Bayesian model
# Equations for the model
model <- alist(
  Sp..Atk ~ dnorm(mu, sigma),   # Likelihood
  mu <- alpha + beta * Attack,  # Linear model
  alpha ~ dnorm(100, 50),         # Intercept Prior
  beta ~ dnorm(1, 0.5),          # Slope Prior
  sigma ~ dunif(0, 200)              # Standard Deviation Prior
)

# Fit the model with Laplace Approximation
fit <- quap(
  model,
  data = data
)

# Summary
precis(fit)
# Graph the model and uncertainty lines
plot(data$Attack, data$Sp..Atk, pch = 16, col = "blue", xlab = "Attack", ylab = "Sp..Atk")
post <- extract.samples(fit)
for (i in 1:100) {
  curve(post$alpha[i] + post$beta[i] * x, add = TRUE, col = alpha("red", 0.2))
}

Respecto a lo que nos da el comando precis: 1. Podemos decir que el intercepto parece razonable considerando el rango entre 0 y 200. La media de 41.31 nos indica que cuando el ataque es 0, el ataque especial esperado es 41.31. Esto es razonable, ya que es esperable que el ataque especial sea mucho mayor que el normal. Además, el intervalo de credibilidad parece razonable para ese mismo valor.

  1. La pendiente parece mostrar un valor bajo entre 0.35 y 0.45. Esto nos indica que la relación entre attack y special attack es positiva, pero attack no es tan influyente en el valor de special attack y hay otras variables que influyen en esta última.

  2. Una alta media para la desviación estándar, lo que nos indica que quizá, el modelo no captura completamente la variación de los datos. Esto sugiere, al igual que con la pendiente, que hay otros factores o variables relevantes que afectan a special attack en los datos.

Lo siguiente es el código que grafica el MAP de la regresión lineal con los intervalos del 95% que se tienen sobre la media y sobre el valor predecido de Sp..Atk. Los resultados son buenos, ya que el intervalo contiene a la mayoría de los valores y las líneas de incertidumbre. El modelo logra ser un buen predictor de los valores estudiados.

ggplot() +
  geom_point(aes(x=poke.df$Attack, y=poke.df$Sp..Atk), color="red") +
  geom_line(aes(x=, y=), color="orange") +
  geom_ribbon(aes(x=, y=, ymin=, ymax=), alpha=0.3, fill="orange") +
  geom_ribbon(aes(x=, y=, ymin=, ymax=), alpha=0.1, fill="orange")

# Generate predictions on the fitted model
attack_seq <- seq(min(data$Attack), max(data$Attack), length.out = 100)

predicted_mean <- link(fit, data = list(Attack = attack_seq))

mean_pred <- apply(predicted_mean, 2, mean)
ci_mean <- apply(predicted_mean, 2, PI, prob = 0.95)

predicted_values <- sim(fit, data = list(Attack = attack_seq))

ci_pred <- apply(predicted_values, 2, PI, prob = 0.95)

plot(data$Attack, data$Sp..Atk, pch = 16, col = "blue", xlab = "Attack", ylab = "Sp..Atk")

# Graph lines and credible intervals
lines(attack_seq, mean_pred, col = "red", lwd = 2)

polygon(c(attack_seq, rev(attack_seq)),
        c(ci_mean[1, ], rev(ci_mean[2, ])),
        col = alpha("red", 0.2), border = NA)

polygon(c(attack_seq, rev(attack_seq)),
        c(ci_pred[1, ], rev(ci_pred[2, ])),
        col = alpha("blue", 0.2), border = NA)

legend("topleft", legend = c("Media predicha (MAP)", "IC 95% (media)", "IC 95% (valores)"),
       col = c("red", alpha("red", 0.2), alpha("blue", 0.2)), lwd = c(2, NA, NA), pch = c(NA, 15, 15),
       pt.cex = 2, bty = "n", cex = 0.8, bg = "white")

 

A work by CC6104